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We study the time evolution of perturbations in spatially extended chaotic systems in the presence of quenched 
disorder. We find that initially random perturbations tend to exponentially localize in space around static pinning 
centers that are selected by the particular configuration of disorder. The spatial structure of typical perturbations, 
8u(x,t), is analyzed in terms of the Hopf-Cole transform, h(x,t) = In \Su(x, t)\. Our analysis shows that 
the associated surface h(x, t) self-organizes into a faceted structure with scale-invariant correlations. Scaling 
analysis of critical roughening exponents reveals that there are three different universality classes for error 
propagation in disordered chaotic systems that correspond to different symmetries of the underlying disorder. 
Our conclusions are based on numerical simulations of disordered lattices of coupled chaotic elements and 
equations for diffusion in random potentials. We propose a phenomenological stochastic field theory that gives 
some insights on the path for a generalization of these results for a broad class of disordered extended systems 
exhibiting space-time chaos. 

PACS numbers: 05.45 Jn, 05.45.Ra, 05.40.-a 



I. INTRODUCTION 

Spatially extended chaotic systems (SECS) are of great im- 
portance for the understanding of fundamental problems in 
deterministic many-particle systems, including hydrodynam- 
ics and turbulence III |2[] or weather forecasting yfl. These 
systems exhibit dynamical instabilities that can be quantified 
by the spectrum of Lyapunov exponents and their correspond- 
ing Lyapunov vectors (LVs) JUSM!,!^]. 

Recently, the application of some tools and concepts bor- 
rowed from the physics of nonequilibirum statistical systems 
has been shown to be a very promising line of research to 
study certain aspects of SECS. In particular, it has been ob- 
served that (after a suitable logarithmic transformation) the 
evolution of an infinitesimal perturbation in deterministic 
SECS can be generically described as a scale-invariant rough 
surface 

JHEI3IMII1- In this surface picture, the erratic fluc- 
tuations of the system, due to the chaotic dynamics, are inter- 
preted as noise. Remarkably, it has been shown Oils [Oh that 
in many cases the associated surfaces belong to the Kardar- 
Parisi-Zhang (KPZ) 11211 universality class of nonequilibrium 
surface roughening. Nonetheless, the probability distribution 
of the randomness generated by the chaotic trajectory can be 
crucial. In this regard, the existence of long-range correla- 
tions |0] or a fat tail of the distribution [ 10] may change the 
universality class observed in the surface growth picture. 

Typically the above discussed studies deal with homoge- 
neous systems, where only identical elements are coupled, of- 
ten in a diffusive manner. However, much less is known about 
the evolution of chaotic perturbations in inhomogeneous ex- 
tended systems, where the many coupled elements are either 
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different or the coupling itself is a quenched random field that 
varies along the system. Such conditions naturally arise in 
some important applications, e. g. in regional weather fore- 
casting, where inhomogeneity is present in the form of ex- 
plicit surface-topography dependent terms JH [3, H, [l5h in 
the dynamic equations. 

Inhomogeneous SECS can indeed demonstrate rather un- 
usual properties, like for instance taming of spatiotemporal 
chaos induced by disorder 01611 . disorder-enhanced synchro- 
nization lfl7i Il8ll . and avoided crossing and level repulsion 
similar to that ocurring for energy eigenvalues in disordered 
quantum systems IU9I1 . 

In this paper, we focus on the dynamics of perturbations 
in disordered SECS, which has been little investigated in the 
literature. We analyze here a rather simple but enlighten- 
ing model consisting of a lattice of coupled chaotic elements 
whose parameters are randomly distributed or the coupling 
among them is a quenched random variable. We show that 
long-range temporal correlations, induced by the quenched 
disorder, can cause perturbations to strongly localize. The 
strong localization observed here is essentially different from 
the dynamical localization that has been previously reported 
for homogeneous SECS ||7lLL!|]. In the latter, the positions of 
the localization centers keep fluctuating in space, while in the 
presence of quenched disorder we find that the positions of 
the localization centers are fixed for a given disorder realiza- 
tion. We analyze typical perturbations, Su(x, t), by making 
use of the Hopf-Cole transform, h(x,t) = In \8u(x, t)\, and 
mapping to the equivalent surface growth problem. Strong lo- 
calization leads to the formation of faceted structures in the 
corresponding surface picture. These Lyapunov surfaces also 
exhibit a coarsening behavior and anomalous kinetic rough- 
ening, akin to some nonequilibrium growing surfaces. The 
connection of the problem of the propagation of errors in dis- 
ordered SECS with the problem of diffusion in random po- 
tentials is also discussed. We propose a phenomenological 
stochastic field theory pointing toward a generalization of our 
results for a broad class of disordered extended systems ex- 
hibiting space-time chaos. 
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H. DESCRIPTION OF THE MODELS 

We focus our numerical study on coupled-map lattices with 
quenched disorder. We consider L coupled chaotic maps 
u x (t), with x = 1, 2, ■ • -L, following the evolution equation 

u x (t + 1) 

&X.X-\-l 

+ (1 - 2a x>x )f x (u x (t)), 

where the a x x ±\ are the disordered nearest-neighbor cou- 
pling constants, a x>x is the on-site contribution, and f x is 
the local nonlinear map at site x. Here we choose f x (g) = 
c x (l/2 — \g mod 1 — 1/2|), which is a periodic continua- 
tion of the tent-map. This map shows chaotic behavior for 
c x G (1, oo), with a Lyapunov exponent that grows logarith- 
mically with c x . This choice allows one to study extended 
coupled systems with a spatial distribution of Lyapunov expo- 
nents by allowing c x to take randomly distributed values along 
the lattice. 

Almost any initial random perturbation Su x (t — 0) will 
grow in magnitude and develop space-time correlations while 
propagating along the system, quickly aligning with the most 
unstable direction in tangent space, the so-called main LV: 



5ux(t + 1) xf x+ x (u x +x (t))8u x+ x (t) 

+o- x , x -xf x -x( u x-i(~ t ))5u x -x(t) 
+(l-2a x , x )f' x (u x (t))5u x (t), 



(1) 



where f x (g) is just the derivative of the local map f x {g) with 
respect to its argument g. 

Numerical integration of the tangent space equations (fl~|i 
looks apparently rather simple, however, a caveat is in or- 
der. We found that the introduction of quenched disorder leads 
to the appearance of very large differences among error field 
8u x (t) values at certain sites of the system. The reason for 
this is the quenched nature of the disorder so that, if a large 
value of the random variable is assigned to a given site, it 
will continue giving high contributions for all times. In fact, 
these differences can become so large during the simulation 
that, if we were to naively integrate (fTJ and just multiply by 
some global factor to avoid numerical overflow, the perturba- 
tion field at those sites where the perturbation values are small 
would soon be considered as zero by the computer due to ac- 
curacy limitations. It is rather simple to overcome this tech- 
nical problem by avoiding calculating the 5u x {t) directly, but 
rather computing the quotients -d x (t + 1) = 5u x (t + l)/Su x (t) 
and ip x {t) = 5u x+ i(t) / 5u x {t) instead. We therefore rewrite 
CD as 

— &X,X-\- 

+a x , x -xf x -x( u x-i(t))/<Px-x(t) 
+ {l-2<j x>x )f' x {u x {i)) 



Px(t + 1) = <p x (t) 



Now we can numerically integrate this pair of equations in- 
stead of Eq. CQ. This useful numerical trick resolves the prob- 
lem of possible overflows during the simulations. 




FIG. 1: The coupling configurations of one site to its neighboring 
sites and viceversa for (a) the random barrier (RB) and (b) random 
trap (RT) model. 



We have studied three different scenarios for a disordered 
couple-map lattice system. In the first scenario (model A) we 
consider the case in which the couplings are identical along 
the system, but coupled elements are inhomogeneous. In the 
simplest setting we can model this situation by drawing the 
local map constant c x from a uniform distribution, so we have 



Model A: 



x.y 



e u(i,2) 

= 1/3 



(2) 



where y G {x — 1, x, x + 1}. Note that we chose a homoge- 
neous democratic coupling (a x ,y — 1/3) all along the system. 
In this situation individual elements are more chaotic or less 
chaotic depending on the corresponding value of c x . Larger 
values of c x lead to a locally faster growth of the perturba- 
tion at site x. Strictly speaking, this is only true for uncoupled 
maps. Things are a much more involved here since the maps 
are coupled. Thus, the behavior of the perturbation at one site 
does not only depend on the map at that specific site, but also 
on the dynamics of its neighborhood. 

A different source of disorder we have explored is the ex- 
istence of a quenched random coupling, while the individual 
maps are all identical. This corresponds to systems with a 
disordered diffusion coefficient. From the classical theory of 
diffusion in disordered media we know that two different sym- 
metries are of interest i20ll2lll . These two models are sketched 
in Fig. [T]and described in detail below. 

On the one hand, we have the random barrier (RB) 
model luOl I21I1 . Physically this corresponds to the existence 
of a random potential barrier between every two neighboring 
sites so that the diffusion coefficient from site x to site x + 1 
is equal to that from site x + 1 to site x (see Fig.[T|). In other 
words, the disorder is associated with the bonds connecting 
neighboring sites on the lattice. We call this configuration 
model B: 



Model B: 





= 2 


&x,x-l 


= £x 


Cx,x+1 


= £x+l 


&x,x 


= (e x -\ 



(3) 



l)/2 



where the disorder is drawn from a uniform distribution, e x € 
U(0, eo), and eo being an arbitrary parameter that we take as 
1/3 unless otherwise stated. 

On the other hand, some disordered systems are better de- 
scribed by the so-called random trap (RT) model J20I |2lll . In 
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FIG. 2: Lyapunov spectra for models A, B and C in a lattice of 
L — 256 coupled tent maps after average over ten different disorder 
realizations. 



this case a static disordered diffusion coefficient is assigned to 
each lattice site x instead to bonds. This corresponds to our 
model C, which is defined by 



Model C: 




2 

£x-l 
(x + l 



(4) 



where the disorder is uniformly distributed, e x 6 U(0, eo), 
and again eo can be any number in the interval (0, 1/2), here 
we use eo = 1/3. 

Typical Lyapunov spectra corresponding to the three disor- 
dered models introduced here are shown in Fig.|2]for a system 
of size L — 256 and averages over 10 different disorder real- 
izations. One can see that, depending on the model, one has 
30-80 % of the spectra in the region A > 0. Note that Lya- 
punov spectra are self-averaging, therefore spectra for a single 
realization differ from the average in small fluctuations, which 
for large enough systems are negligible. 

Some remarks are now in order. First, let us stress that these 
two disordered diffusion configurations, namely RB and RT, 
were not chosen arbitrarily, but are in fact all existing phys- 
ically meaningful ways to introduce disorder in the diffusive 
couplings [20, 21]. Second, it is important to remark that in 
this paper we focus on systems in the presence of weak disor- 
der, where the probability density of having zero disorder at 
any given site, V(e = 0), is bounded. We do leave out of our 
study the case of strong disorder, V(e — > 0) ~ e~\"\, where a 
fraction of the system sites can effectively act as sinks, which 
is known to have a great impact on the asymptotic transport 
properties I20ll2lll and will be discussed elsewhere. 

In principle, classical results on the problem of diffusion 
in random media 1I20I I2H1 may also be invoked to argue that 
the details of the disorder distribution should be irrelevant in 
the sense that weak disorder gives rise to the same asymptotic 
dynamical behavior independently of the detailed form of the 
disorder distribution. The same conclusion applies in the case 
of strong disorder, where the asymptotic dynamics is fully de- 
termined by the form of the divergence at the origin, i.e. the 
value of the exponent v. In the following we focus on the case 
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FIG. 3: Typical evolution a Lypunov vector surface h(x, t). One can 
observe the emergence of facets that grow in time until, in the long 
time limit, the whole interface is spanned by one single facet. Tiny 
fluctuations can be observed at small scales. The data correspond to 
model C in a lattice of L = 1024 maps. 



of weak disorder and a uniform distribution V(e) = U(0, eo) 
with the parameter eo = 1/3 is used to exemplify our results. 
However, in the presence of random multiplicative terms the 
critical behavior might be affected by the detailed form of the 
disorder distribution 122]. 



III. SURFACE EVOLUTION AND COARSENING 

As occurs in the case of homogeneous systems, the dynam- 
ics of the main LV can be conveniently described in terms 
of the equivalent surface picture |0> H B OH EH] • The rea- 
son is that the logarithm of the perturbation field turns out 
to be scale-invariant, so that correlations have the form of 
power-law functions and critical exponents can then be used to 
characterize the space-time structure of perturbations. There- 
fore, in this paper we are interested in the scaling proper- 
ties of the rough surface defined by the Hopf-Cole transform 
of the perturbation h(x,t) = \n\5u x (t)\ = \n\du x (0)\ + 

Et=iM<WI- 

Space-time scaling properties of chaotic perturbations in 
inhomogeneous systems turn out to be very different from 
those in homogeneous systems, even in the case of weak dis- 
order studied here. As can be immediately seen in Fig. [3] 
a first observation is the patterned structure of the LV sur- 
face, which is visible to the naked eye and contrasts with 
the KPZ-type morphologies observed in homogeneous sys- 
tems 10, EI]- This already indicates that significant differences 
are to be expected in the scaling properties of inhomogeneous 
systems. Even weak disorder induces self-organization in a 
triangular structure, which reflects the strong spatial local- 
ization of the perturbation in the form of a exponential pro- 
file around some strong pinning centers. Note that a facet, 
h(x,t) = h(xo) — s(t)\x — xq\, with a cusp at x = xo and 
slope s(t) corresponds to an exponential profile of the per- 
turbation, Su x (t) oc exp[— s(t)\x — xq\], around the pinning 
center xq. 

In Fig. |3] we plot a typical evolution of the main LV ac- 
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FIG. 4: (Color online) The evolution of the lateral mound size p(t) 
is plotted for the three models. The straight lines are plotted to guide 
the eye and have slopes 0.78 (dotted), 2/3 (dashed) and 1/2 (solid), 
respectively. The curves for models A and B have been shifted for 
better visibility. 
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FIG. 5: (Color online) The spectral power functions of the 
surfaces obtained for model A are plotted for times t = 
2 12 (a),2 14 (b),2 16 (c),2 18 (d). The straight line corresponds to 
a s = 1.5. In the inset we show the data collapse according to the 
scaling ansatz ((6} and l[7}, using a — 1.03 and z — 1.3. The dashed 
and the solid straight line have slopes -1 and 3, respectively. The 
plotted data correspond to model A for a system of size L = 16384 
averaged over 100 disorder realizations. 



cording to Eq. ([l} in a system of size L — 1024 for a given 
disorder realization in the case of model C. We observe that 
at short times, the surface h(x, t) is composed of triangular 
facets of varying sizes which grow in time until, in the long 
time limit, the whole interface is formed by just one single 
facet. Therefore, as time evolves the size of the triangles in- 
creases and number of triangles diminishes. This surface dy- 
namics corresponds to non-equilibrium 'coarsening'. Similar 
triangular morphologies are obtained for the three models A, 
B and C introduced above. Close inspection and analysis of 
critical exponents shows that the three models actually belong 
to different universality classes. 

In order to characterize the growth of the faceted pattern 
we have measured the lateral mound size p(t), which gives 
the coarsening length or typical length scale of the instabil- 
ity. This is usually done by calculating the slope-slope cor- 
relation function (Vh(x, t)Vh(x + r, t)) and measuring the 
distance, p(t), at which this correlation crosses zero. We cal- 
culate the surface gradient as the centered discrete derivative 
Vh(x,t) = [h(x + l,t) - h(x - l,t)]/2. In Fig.g]the lat- 
eral mound size, p(t), for the three models is shown. After 
a short initial period, pit) grows as a power law in time with 
exponents 9 B = 0.67 ± 0.03 and 9 C = 0.48 ± 0.03, for mod- 
els B and C, respectively. At long times p(t) saturates, as the 
surface becomes dominated by a single facet and the coarsen- 
ing length becomes comparable to the system size. For model 
A the intermediate region does not seems to be following a 
simple power-law. Nonetheless, we have fitted an exponent 
6a = 0.78 ± 0.05 for the sake of comparison. The different 
scaling exponents for the coarsening length already indicate 
that the three models possibly belong to different universality 
classes. 



IV. SURFACE ROUGHNESS SCALING 

Faceted surfaces exhibit anomalous scaling properties in 
the sense that local and global fluctuations may scale with dif- 
ferent scaling functions (and exponents), which can be cast 
in the generic dynamic scaling form introduced by Ramasco 
et. al. |23[|. In this section we briefly describe the generic 
scaling theory for surface kinetic roughening that will be used 
later on to analyze our numerical data. 

The roughness of scale-invariant surfaces corresponds to 
the fluctuations of the surface height. These fluctuations can 
be computed either locally or globally as follows. On the 
one hand, the global roughness exponent a can be obtained 
from the scaling behavior of the global width W(L, t) = 

([h(x, t) — /i(t)] 2 ) 1 / 2 . Here, the overline denotes an average 
over all sites x in a system of size L and brackets denote the 
average over different realizations. For scale-invariant sur- 
faces one expects the global width to scale as W(L,t) = 
t a l z Q{Llt x l z ), where Q{u) is a scaling function that becomes 
constant for u 3> 1 and decays as ~ u a for u <C 1. The 
roughness exponent a and dynamic exponent z characterize 
the scaling behavior of the global surface fluctuations. 

On the other hand, one can measure the local roughness ex- 
ponent aioc, which is defined via the scaling behavior of the 
local width w(l, t) = {([h{x, t) - (h)^)} 2 )^ 1 / 2 , where (. . 
denotes an average over a; in a window of size I, The local 
width scales as w (l,t) - t a / z g A (l/t^ z ), where the scaling 
function Qa(u) has a similar asymptotic behavior as Q(u), 
but with a local anomalous exponent, ~ u Qloc for u <C 1. 
So that in the stationary regime (for t 3> L z ) one obtains 
Wst&t(l,L) ~ / Q i°c£a-aioc i n (jj e cases where local and 
global roughness exponents do not coincide the scaling is said 
to be anomalous. 

Following Ramasco et. al. j23ll . in order to correctly clas- 
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sify the different forms that scaling can take it is convenient 
to introduce a third roughness exponent, namely the so called 
spectral roughness exponent a s . This is defined in terms of 
the structure factor (or power spectrum in k space), 



S(k,t) = (h(k,t)h(-k,t)), 



(5) 



where h(k, t) is the spatial Fourier transform of h(x, t) — h(t). 
Ramasco et. al. [23] showed that scale-invariant roughening 
in d + 1 dimensions is fully described in the following scaling 
ansatz: 



S(k,t) = k'^+^sikt 1 ^), 



(6) 



z being the dynamical exponent that connects temporal- scales 
r and length-scales I according to t ~ l z , and the nontrivial 
scaling function takes the form 



s{u) 



u 2a+d if u < 1. 



(7) 



Here a is the global roughness exponent. With the help of 
these three exponents one can distinguish the four different 
possible types of scaling that any scale-invariant rough surface 
can exhibit, 



if a s < 1 => a loc = a s 
if a s > 1 =>• ai oc = 1 




Family- Vicsek 
intrinsic 
super-rough 
faceted. 



(8) 

Note that only in the case where standard Family-Vicsek scal- 
ing is valid, which means a = ai oc — a s , the surface is 
self-affine. In the other cases this will in general not be true, 
i. e. local and global fluctuations will still exhibit dynamical 
scaling, but will do so with different exponents. 

Although Family-Vicsek scaling is the best known scaling 
class i24ll . over the last ten years experimental studies in a va- 
riety of systems including growth of thin-films, electrodepo- 
sition, fracture or fluid imbibition II25[ bd |271 l28l 29 l 30l bit] 
as well as theoretical studies |23l l^mflMM IMPHM 
M m S3 m B have confirmed the existence of intrinsic 
anomalous and super-rough kinetic roughening, as well as the 
validity of the generic scaling theory. 

We now use this general scaling ansatz to study the scal- 
ing properties of our system. In order to do so, we have cal- 
culated the structure factor corresponding to the Hopf-Cole 
transforms of the perturbations at various times. Figures [5] 
[6] and [7] show the structure factor for models A, B and C, re- 
spectively. We observe a power-law decay as ~ fc-( 2a =+ 1 ) 
for momenta limited by a large scale and a small scale cut- 
off. At large length scales the power-law behavior is limited 
by the correlation length £(t) ~ i 1 / 2 , as expected from dy- 
namical scaling. The cutoff at short lengths also grows in 
time, suggesting the existence of a different dynamics at small 
scales, which will be studied in detail in Sec. V. Note that the 
curves S(k, t) shift downward for increasing times, indicating 
the presence of anomalous scaling with a s > a. Specifically, 
one expects S(k, t) ~ k 2aa+1 t( a ~ aiac ^ z . We have measured 




FIG. 6: (Color online) The spectral power functions of in- 
terfaces obtained for model B are plotted for times t = 
2 12 (a) , 2 14 (b) , 2 16 (c) , 2 18 (d) , 2 20 (e) . The straight line corresponds 
to q s = 1.5. In the inset we show the data collapse according to scal- 
ing ansatz ((6} and 10, using a — 1.0 and z — 1.5. The dashed and 
the solid straight line have slopes -1 and 3, respectively. The plot- 
ted data corresponds to model B for a system of size L — 16384 
averaged over 100 realizations. 



a s = 1.43 ± 0.05, 1.47 ± 0.05 and 1.45 ± 0.05 for mod- 
els A, B and C, respectively. This should be compared with 
the value of a s = 3/2 that can be obtained analytically for 
interfaces consisting of random smooth facets 12311 . In the in- 
sets of Figs.|5]|6]and|7]we show a data collapse according to 
the generic scaling ansatz in Eqs. © and (0, obtaining val- 
ues for the roughness and dynamical exponents a and z. As 
can be verified, all three models show scaling properties well 
described by the faceted interface scaling class in Ramasco's 
classification, Eq. ([8]), with a roughness exponent awl for 
all the three models. Nonetheless, all three systems yield dif- 
ferent values for the dynamical exponent, z\ = 1.3 ± 0.1, 
z B = 1.50 ± 0.03 and z c = 2.00 ± 0.03, again indicating 
that the three models belong to different universality classes. 
These values of z define the typical length scale that correla- 
tions have spread along the system up to time t, ~ t 1 ^, 
and are in good agreement with the inverse of the exponents 
(8 = 1/z) that we obtained from the analysis of the coarsen- 
ing length discussed in Sec. Ill for all the three models studied 
in this paper. 



V. SEPARATING FACETS FROM FLUCTUATIONS 

The scaling analysis of the surface fluctuations presented 
in the preceding sections has shown that the LV surface self- 
organizes in a characteristic triangular pattern for all the three 
models of disorder. A closer inspection, in particular the com- 
putation of the structure factor, has also revealed the exis- 
tence of a cutoff at small scales such that length scales be- 
low that point obey a different dynamics. This suggests the 
existence of another dynamical process taking place at small 
scales which is different from the mechanism responsible for 
the large scale faceted structure. 

In this section we show that one can actually separate two 
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FIG. 7: (Color online) The spectral power functions of in- 
terfaces obtained for model C are plotted for times t — 
2 12 (a), 2 14 (b), 2 16 (c), 2 18 (d), 2 20 (e). The straight line corresponds 
to a a — 1.5. In the inset we show the data collapse according to scal- 
ing ansatz $6^ and Q, using a = 1.0 and z — 2.0. The dashed and 
the solid straight line have slopes -1 and 3, respectively. The plot- 
ted data corresponds to model C for a system of size L — 16384 
averaged over 100 realizations. 
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FIG. 8: For a segment of 1024 maps on a system of L = 16364 
we exemplify the separation of the surface h(x) into the triangular 
pattern h°(x) and the difference y(x) — h(x) — h°(x), which gives 
the small scale fluctuations. 



different contributions to the LV surface height: The large 
scale pattern, h°(x, t), formed by facets of constant slope, and 
a randomly fluctuating term, y(x, i), which becomes the dom- 
inant one at small scales. Therefore the surface profile at any 
given time can be expressed as the sum of the two independent 
components as h(x,t) = h°(x,t) + y(x,t). This separation 
can be carried out numerically as follows. For a given time 
the interface h(x, t) is smoothed to remove local maxima and 
minima corresponding to local fluctuations upon the triangu- 
lar structure. This is done by replacing the heights at every 
site by a spatial average over some arbitrary region around 
each site. Then, the locations of the cusp sites and valley sites 
of the smoothed curve are easily identified and the faceted 
pattern h°(x) is defined as the set of straight lines connecting 
the cusp points to the neighboring valley points. The interface 
corresponding to the fluctuations is then obtained by taking 
the difference y(x,t) = h(x,t) — h°(x,t) as exemplified in 
Fig.E 



FIG. 9: The power spectra of the interface h(x) (light gray dots), 
the triangular structure h°(x) (dark gray lines), and the difference 
y(x) — h(x) — h°(x) (gray lines) are plotted for two different times 
t& < tt>- The straight solid line corresponds to a" = 1.5, while the 
straight dashed line corresponds to af = 0.45. 

Interestingly, one can observe that the triangular pat- 
tern, h°, and the random fluctuation components, y(x,t), 
are uncorrected. This can be shown by checking that 

(h°(k, t)y(— k, t)) = 0, which is equivalent to proving the 
identity 

S(k, t) = (hP(k, t)hP(-k, t)} + (y(k, t)y(-k, <)). (9) 

In fact, the numerical data presented in Fig. [9] confirm that 
both components are actually uncorrected. Figure [9] shows 
the structure factor of the faceted pattern h°(x,t), the local 
fluctuations y(x, t), and the complete interface h(x, t) at two 
different times. In particular, the data shown correspond to 
model B, but we obtained identical results for for the other 
two models. 

It becomes apparent that the anomalous scaling of the sur- 
face stems from that of the faceted structure. Since the vertical 
scale in Fig. [9] is logarithmic, the sum in Eq. (O will essen- 
tially be dominated by the largest of the two terms, as one can 
easily see in Fig. [9] Therefore, the large scale behavior of the 
complete interface is totally dominated by the scaling behav- 
ior of the pattern component. On the other hand, the structure 
factor S(k, t) should crossover to that of the local fluctuations, 
(y(k,t)y(—k,t)), at short wavelengths. The structure factor 
of the pattern component shifts downwards with time because 
a < a s . This is not the case for the structure factor of the local 
fluctuations, which scale according to Family-Vicsek's, where 
a = a s « 0.45 is fulfilled. Thus the length scale at which the 
local fluctuations dominate the structure factor should grow 
with time. Whether the asymptotic scaling behavior of the 
facets can be observed in the limit t — > oo for a finite system 
depends on the relative strength of the local fluctuations. 

VI. LYAPUNOV VECTORS 

So far we have studied the dynamics of random infinitesi- 
mal perturbations, which is equivalent to study the dynamics 
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of the most unstable direction in tangent space, i. e. the first 
LV. We now devote this section to briefly describe the dynam- 
ics of further unstable directions that growth at slower rates. 

Let us start by introducing a few key quantities that shall 
be used in our analysis. The evolution equation for an ini- 
tially random infinitesimal perturbation, Eq. (|T), can be writ- 
ten in vectorial form as Su(t + 1) = J[u(t)]5u(t), where 
J[u(t)] is just the jacobian evaluated on the trajectory u(t) 
at time t. Since infinitesimal perturbations evolve linearly 
there exists a linear operator M(<2,ii) such that Sufo) — 
M(t a ,ti)«u(*i). given by M(t 2 ,h) = UuZl^HU)]- Os- 
eledec's theorem [44] implies that the (symmetric) limit oper- 
ator limi^ooIM^^OM*^,^)] 1 / 2 ^ 2- * 1 ), where M* de- 
notes the adjoint of M, does exist and the logarithms of its 
eigenvalues are the Lyapunov exponents A n (n = 1, 2, ■ ■ • , L) 
and the eigenvectors are the so-called backward Lyapunov 
vectors b„ [0, 0. In what follows, we will consider stan- 
dard ordering of the Lyapunov exponents, Ai > A2 > . . . > 
. . . Xl, and the corresponding vectors. We use standard nu- 
merical techniques to calculate Lyapunov exponents and vec- 
tors HHH. 

As we have discussed in Sees. Ill and IV, in inhomogeneous 
systems the main LV, the one corresponding to the most un- 
stable direction, converges to a faceted structure asymptoti- 
cally dominated by one single facet in the limit t — > 00. In 
Fig. [TO] we plot the LV surfaces h n (x,t) — ln(\b n (x,t)\), 
n G {1,2,3,4,5}, corresponding to the first five backward 
LVs for model B at some arbitrary time t during evolution in 
a system of 512 coupled maps. The surface profiles in this 
plot are arbitrarily shifted in the vertical direction to aid com- 
parison. In the long time limit the interface hi(x, t), the first 
LV surface, naturally coincides with the interface h(x, t) cor- 
responding to an arbitrary initially random perturbation. It 
should be observed that not only the first LV, but each surface 
h n (x, t) corresponding to the n-th LV indeed converges to a 
structure dominated by a single triangle as well. This is actu- 
ally the case for all the inhomogeneous models studied in this 
paper. 

Interestingly, the location of the cusp of the triangles, where 
the LVs spatially localize, is not arbitrary. It is interesting to 
compare the positions of the global maxima of the interface 
profiles with the values of the quenched disorder at those sites. 
For model C, i. e. the RT case, the site where the first LV 
takes its maximum value coincides with the site x at which the 
coupling e x takes its lowest value, the site where the second 
LV takes its maximum coincides with the site where e x takes 
its second lowest value, and so forth. This is easy to interpret: 
since errors get trapped at sites with lower diffusion, the result 
is that the LVs strongly localize and utterly get frozen at those 
sites. 

A very similar correlation effect between error localization 
and disorder can be seen for model B, e. the RB case. Note 
that in the case of RB disorder the diffusion coefficient at a 
given site x is the sum e x +i + e x of the local diffusive cou- 
plings connecting x to any of the two neighboring sites. In 
fact, we find that the site where the first LV takes its maxi- 
mum coincides with the site x where the diffusivity e x+ i + e x 
takes its lowest value, the site where the second LV takes its 




X 



FIG. 10: The surfaces h n {x) corresponding to the first five 
Lyapunov-vectors are plotted (from bottom n = 1 to top n = 5). 
Note that all the interfaces converge to a triangular structure. The 
data presented were obtained for model B and have been shifted for 
better visibility. 



maximum value coincides with the site where e x +i + e x takes 
its second lowest value and so forth. Thus, the on-site growth 
velocity is again controlled by the probability for the pertur- 
bation to get trapped at that site. 

One would expect that a similar association of LVs maxima 
with disorder extrema should also exist for model A. How- 
ever, things are not that straightforward in this case. We recall 
that in model A the disorder is introduced in the map con- 
stants and not in the diffusive couplings, which are identical 
along the system. We think that due to the dominance of small 
scale fluctuations up to quite large length scales for model A 
(see Fig. [5]l as compared with the situation for models B and 
C, the maxima get smoothed out and an extended neighbor- 
hood, rather than just one single site, controls the final position 
where the LVs get frozen at long times. In other words, local 
synchronization of neighboring sites over short scales makes 
it difficult to simply identify the region where, at a coarse- 
grained level, the vectors will get asymptotically trapped in 
the case of disorder of type A. 

VII. TORWARD A PHENOMENOLOGICAL STOCHASTIC 
FIELD THEORY 

The evolution of perturbations in homogeneous SECSs is 
known to be described in statistical terms by the multiplicative 
stochastic equation 

d 

—Su = d xx 5u + £(x,t)6u, (10) 
at 

where £(x,t) is a noise term that accounts for the chaotic 
fluctuations along the trajectory and is assumed to be delta- 
correlated in space and time, (£(x, t) £[&', £')) = 2a8(x — 
x') 5(t — t'). Pikovsky and Kurths (|7[] and Pikovsky and 
Politi [8] proposed this equation as the proper candidate for 
modeling the dynamics of the first Lyapunov vector from a 
statistical perspective. They showed that it actually repro- 
duces the main statistical properties of SECS in a variety of 
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systems [8]. This equation mimics the linear equations in tan- 
gent space for the dynamics of infinitesimal perturbations in 
spatio-temporal chaotic systems. 

An important feature of Eq. ( TTOb is that, under a Hopf-Cole 
transformation, h = In \Su\, it maps into the KPZ equation: 

|U = (d x h) 2 + d xx h + t{x,t), (ID 

at 

which ultimately justifies why the log-transformed (main) 
Lyapunov vector of many spatiotemporal chaotic systems 
scales in space and time as a KPZ surface 10, HI |H [Tfl [Till . 
Interestingly, Eq. (TTOb also appears in the context of the classi- 
cal problem of the directed polymer in a random potential fiTtl 
(also see Ref. [48] for an detailed review). 

We now explore extensions of this type of Langevin equa- 
tions to the case of disordered systems. A phenomenological 
stochastic field theory of this kind should be helpful in order to 
identify the essential symmetries and mechanisms that deter- 
mine the dynamics of perturbations in chaotic disordered sys- 
tems. The existence of scale invariance strongly encourages 
this approach. Therefore, microscopic details of the models 
are expected to be irrelevant at a sufficiently coarse-grained 
scale, as occurs in the case of homogeneous systems. 

In the case of Model A, Eq. (fTJ can be considered as a 
straightforward discrete version of 

d 

— 5u = d xx [d u f(x,u)Su] + [d u f{x,u) - l]Su. 
at 

It is important to note that the noisy term d u f(x,u) = 
d v f(x,y)\ u ( x t ) depends implicitly on time through the tra- 
jectory u(x,t) and so, it includes both an annealed and a 
quenched component. The origin of the latter is the inhomo- 
geneous character of the coupled elements: Sites are more 
or less chaotic depending on their position x in the system. 
For instance, in the particular case of the logistic map we 
have d u f(x,u) = c(x)[l — 2c(x)u], where c(x) is the local 
map parameter. We argue that in the hydrodynamic limit the 
quenched part of the noise dominates over the annealed one. 
The reason being that the quenched character of the disorder 
generates long-range temporal correlations at each site, which 
are expected to control the dynamics over the short-range cor- 
relations of the annealed random terms. 

These arguments suggest that for type A models the long 
time limit of main LV should be generically described by a 
multiplicative Langevin equation like 

d 

— 5u = d xx [(((x)Su\ + d xx 6u + ((x)6u, (12) 
at 

where ((x) represents the quenched disorder. This equation 
can be further simplified since the disordered diffusion term 
d xx [(C( x )^ u ] is irrelevant as compared with d xx 8u. To be 
precise, in the case of weak disorder that we study here the 
inhomogeneous diffusion term simply renormalizes to an ef- 
fective constant diffusion term [ 60] . Finally we arrive at our 
stochastic field theory proposal for model A type of systems: 



where D is an effective diffusion constant and Q(x) is a 
quenched columnar disorder with delta correlations. 

Interestingly, Eq. ( PT31 appears in the context of diffusion 
in random trappin g/ amplifying quenched potentials 1491 l50l 
I55L 15611 . After the Hopf-Cole transform of the 
perturbation field, 8u — exp([ih/D), one obtains 

d t h{x,t) = d xx h + f i(d x h) 2 + ((x), (14) 

which is a close relative of the KPZ equation, but here the 
additive noise term ((x) is quenched and delta correlated, 
(C( X )C( X ')) — 205(x — x'). This equation has attracted some 
interest in the past, as it also describes the free energy of di- 
rected and undirected polymers of length t in the presence of 
columnar disorder Hi H3 HI HI . 

Only very recently [59], it has been realized that the sur- 
face described by Eq. ( fT4b exhibits anomalous scaling ex- 
ponents induced by the self-organization in a faceted struc- 
ture. In Ref. 15911 we studied Eq. ( fT4b by means of exten- 
sive simulations and determined the critical exponents to be 
z = 1.35 ± 0.05, a = 1.05 ± 0.05 and a s = 1.50 ± 0.05. 
These values are in excellent agreement with the values re- 
ported here for the critical exponents of the main LV surface 
in SECSs of type A (see Fig. 0. This strongly supports our 
claim that, in fact, Eq. ( fT3l is a minimal model that captures 
the essential ingredients dominating the evolution of errors in 
inhomogeneous models of type A 16111 . 

As it has already been discussed in Ref. [59], the introduc- 
tion of an additional annealed noise term in Eq. ( fT4b is ir- 
relevant in the renormalization group sense, but leads to an 
increase of the typical size over which the local random fluc- 
tuations (discussed in Sec. V) dominate. Hence, for a finite 
size system the asymptotic scaling behavior may not be ob- 
servable if the thermal noise amplitude is too large. This is 
of some importance for the situations we are interested in 
because, when considering the evolution of perturbations in 
SECS, there will always be some temporally fluctuating noise 
as explained above. Thus, even if quenched disorder is present 
in the studied system, it may be very weak in comparison with 
the temporally fluctuating one. In that case the behavior de- 
scribed here would only be observable in very large systems. 

Moreover, it should be noticed that, although the roughness 
exponents a, a\ oc and a s seem to be the same within error 
bars for the three models of disorder studied in this paper, the 
dynamic exponent z does take different values. This suggests 
that the disorder symmetries that distinguish models A, B and 
C are indeed relevant and lead to different nonequilibrium uni- 
versality classes. One can try to formulate a stochastic-field 
evolution equation similar to Eq. (fT2l for models B and C, but 
where the disorder enters in a different form into the diffu- 
sion and the multiplicative term in order to reflect the differ- 
ent symmetries of the coupling in these two models. Unfor- 
tunately, we have not been able to find the correct equations 
reproducing the scaling behavior of the LV surfaces for disor- 
dered systems in the class B and C. Similar arguments might 
be used to write down the stochastic equation 



dtSu = Dd xx 5u + £(x)8u, 



(13) 



d t Su = d x [((x)d x 5u] + £(x,t)6u, (15) 
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for model B and 

d t Su = d xx [((x)6u] + £(x, t)5u, (16) 

for model C, where £(x, t) is a noise term. We have inves- 
tigated the scaling properties of Eqs. (T5[ and ( TToT l by means 
of numerical simulations and found that, while these simple 
stochastic models generically lead to faceted surfaces with the 
expected roughness exponents, they are unable to reproduce 
the correct dynamic exponent of disordered SECS in the class 
B nor C. We believe that the reason behind this disagreement 
lies in the existence of correlations between the quenched and 
annealed fluctuating terms that are not incorporated in this 
naive approach. These non-trivial correlations might also ex- 
plain the, at least at first sight, counter intuitive result that dif- 
ferent dynamic exponents are observed in SECS of class B 
and C. Some remarks regarding this point are now in order. 

Classical diffusion in random media involves the study of 
equations like (fT3T > and ( fT6b , but where the multiplicative term 
t)8u is absent. It is well-known [21] that, in the long time 
limit, dynamics of RB and RT models is identical, even for 
strong disorder. Despite the different symmetries involved, 
the large scale physics is the same. However, our numerical 
results indicate that model B (RB symmetry) and model C (RT 
symmetry) in SECS seemingly belong to different universal- 
ity classes. We claim that the multiplicative noise, which is 
generically coupled to the quenched disorder in the systems 
we have studied, breaks the duality between RB and RT type 
of disorder in this case, leading to distinct dynamical expo- 
nents for the two different diffusion configurations. 

VIII. CONCLUSIONS 

We have studied the evolution of infinitesimal perturba- 
tions in inhomogeneous spatially extended systems exhibit- 
ing space-time chaos. Inhomogeneity is introduced by means 
of a quenched disorder. We have considered one-dimensional 
coupled-map lattices as a simple and computationally conve- 
nient model system to analyze some aspects of chaos in the 
presence of disorder. In this regard, three different classes of 
models have been investigated by means of extensive numer- 
ical simulations. In all cases we find a strong localization be- 
havior characterized by an exponential spatial profile around 
some localization centers. In the long time limit the perturba- 
tion concentrates around one single final attracting center. We 



have also studied the second, third, and so on backward LVs 
corresponding to most rapidly expanding directions in tangent 
space and found a correspondence between the localization 
centers and the positions in space corresponding to increasing 
minima of the diffusion in the particular disorder realization. 

Note that the strong localization behavior described here is 
essentially different from the dynamic localization observed 
in the case of homogeneous SECS. In the latter, perturbations 
do also localize on just a few sites, but the position of these 
sites keeps fluctuating in time. However, in the presence of 
quenched disorder the sites where perturbations localize are 
fixed by the corresponding realization of the disorder. 

Moreover, by a standard mapping (Hopf-Cole transform) 
of the perturbation into a growing surface we found that the 
LV associated surfaces self-organize in a faceted structure, at 
variance with what occurs in homogeneous (non disordered) 
systems where one generically obtains a surface in the uni- 
versality class of KPZ. Interestingly, this faceted surface was 
found to exhibit coarsening and anomalous kinetic roughen- 
ing in agreement with previous theoretical predictions Il23ll for 
this type of scale-invariant structures. 

The evolution of infinitesimal perturbations in spatially ex- 
tended chaotic systems with quenched disorder may be de- 
scribed at a coarse-grained level as a diffusion process in a 
random potential. Note, however, that the value of the dy- 
namical exponent z depends on how disorder is introduced. 
The reason for this is possibly the influence of the disorder 
distribution on the temporal dependence of the coarsening, 
as described in Ref. [22] in the context of diffusion (see also 
Ref. B59I1 for a discussion relevant to the present work) 

We think that strong localization and anomalous scaling 
should play a role, e. g. in realistic weather models, where 
quenched disorder is included in the form of inhomogeneous 
boundary conditions representing certain geographical and to- 
pographical conditions. 
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